package Market;

import cern.colt.matrix.*;
import cern.colt.matrix.impl.*;

public class OptimalPowerFlow {
	
	  // Indices for Node Data parameters: form of {number of buses, Penalty coefficient}												 
	  private static final int NInd =       0;  				
	  private static final int PenCoefInd = 1;
	  
	  
	  // Indices for Branch Data parameters: form of {From, To, Branch Capacity, Reactance}
	  private static final int FromInd =  0;
	  private static final int ToInd =    1;
	  private static final int BCapInd =  2;
	  private static final int ReactInd = 3;
	  
	  
	  // Indices for Supply Offer parameters: form of {a,b,CapL,CapU}
	  private static final int aInd     = 0;
	  private static final int bInd     = 1;
	  private static final int CapLInd  = 2;
	  private static final int CapUInd  = 3;
	  														/** number of Nodes */
	  private int K;            
	  														/** number of Branches */
	  private int N;              
	  														/** number of GenCos */
	  private int I;           
	  														/** number of LSEs */
	  private int J;          
	  
	  														/**  Node Data Array: {number of buses, Penalty coefficient} */
	  private double[][] NodeDataA;        
	  														/** Branch Data Array: {From, To, Branch Capacity, Reactance} */	  
	  private double[][] BranchDataA;
	  
	  double[][] BranchDataA2;
	  														/** Hourly loadProfile Array = fixed demand : { each LSE demand over current hour } */	  
	  private double[]   LoadProfileA;     
	  														/** True Supply Offer Array: {a,b,CapL,CapU} */  
	  private double[][] TSupplyOfferA;     
	  														/** Array of Generators' locations {gen1 at node 1, gen2 at node3,.., geni at nodej} */
	  private int[] GenAtNodeA;            
	  														/** Array of LSE locations {LSE1 at node 2, LSE2 at node5,.., LSEi at nodej} */
	  private int[] LSEAtNodeA;            
	  														/** N (1x1) number of branches */
	  private int nBranches;                    
	  														/** K (1x1) number of nodes */
	  private int nNodes;                    
	  														/** Matrix to hold Branch Data:  {From, To, Branch Capacity, Reactance} */
	  private DoubleMatrix2D BranchDataM; 
	  														/** Matrix to hold Node Data: {number of buses, Penalty coefficient}  */
	  private DoubleMatrix2D NodeDataM;  
	  
	  
	  // Input for QuadProgJ (G,a,Ceq,beq,Ciq,biq)
	  														/** Matrix G =blockDiag [U  Wrr] of (I + K − 1)× (I + K − 1) dimension <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private DoubleMatrix2D GM;
	  														/** Vector a=[A1 ... Ai 0 ... 0] of (I + K − 1) dim, where Ai are a's from trues supply offer/pu <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private DoubleMatrix1D aM;
	  														/** The equality constraint matrix Ceq =[ II −Br ] of K × (I + K − 1) dim  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20)*/ 
	  private DoubleMatrix2D CeqM;
	  														/** The associated equality constraint vector  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */ 
	  private DoubleMatrix1D beqM;
	  														/** The inequality constraint constraint matrix <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 20)*/ 
	  private DoubleMatrix2D CiqM;
	  														/** The associated inequality constraint vector <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 21)*/ 
	  private DoubleMatrix1D biqM;
	 
	  private QuadProgJ qpj;
	  
	  // Solution from QuadProgJ
	  private double[] Commitment; // power production quantity
	  private double[] VoltAngleR;  // voltage angle in radians
	  private double[] VoltAngleD;  // voltage angle in degrees
	  private double[] LMP;  // lmp (locational marginal prices) == eqMultiplier
	  private double[] IneqMultiplier;
	  private String [] ineqMultiplierName;
	  private double minTVC; // minTVC = SUM (Ai*PGi + Bi*PGi^2)
	  private double[] BranchFlow; // branchFlow_km = (1/x_km)(delta_k - delta_m)
	  private double[][] bi; // branch index
	  private double sumSquaredAngleDifference;  // SUM(delta_k - delta_m)^2
	  
	  
	  // Intermediate input to form (G,a,Ceq,beq,Ciq,biq)
	  														/** Vector B=[B1 ... Bi] of Ix1 dim, where Bi are b's from trues supply offer/pu */
	  private DoubleMatrix1D B;  
	  														/** An attribute matrix U = diag[2B1 , 2B2 , ··· , 2Bi] of Ix1 dim <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17)*/
	  private DoubleMatrix2D UM;    
	  														/** Reduced weight matrix  <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17)*/
	  private DoubleMatrix2D WrrM; 
	  														/** Vector A=[A1 ... Ai] of Ix1 dim, where Ai are a's from trues supply offer/pu <br>  @see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
	  private DoubleMatrix1D A;  
	  														/** II Matrix represents K Nodes(rows)  and I Generators(columns) located on them */
	  private DoubleMatrix2D IIM;
	   														/** Nx1 matrix  of branches' capacities  */
	  private DoubleMatrix1D BCapM;   
	  														/** Ix1  matrix of generator's true lower capacity */
	  private DoubleMatrix1D CapTLM;  
	  														/** Ix1 matrix of generator's true upper capacity */
	  private DoubleMatrix1D CapTUM; 
	  														/** (N×I) zero O matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */ 
	  private DoubleMatrix2D OtM;
	  														/** (I x K-1) zero O matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21)  */
	  private DoubleMatrix2D OpM;
	  														/** I×I Identity Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	  private DoubleMatrix2D IpM;
	  														/** (NxN)Diagonal Admittance Matrix, where Bkm=1/Xkm <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	  private DoubleMatrix2D DAM;
	  														/** (N x K-1) Reduced Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19)  */
	  private DoubleMatrix2D RAM;
															/** True Supply Offer Matrix: {a,b,CapL,CapU} */
	  private DoubleMatrix2D TSupplyOfferMatrix; 														
	  														/** (KxK) Forms Voltage Angle Difference Weight Matrix <br>@see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */
	  private double[][] VADWM;               
	  														/** (K-1)x(K-1) Reduced Voltage Angle Difference Weight Matrix <br>@see  DC-OPF.IEEEPES2007.JSLT.pdf (pg 17) */
	  private double[][] RVADWM;             
	  														/** (Nx2) Matrix which holds the data for Branches' connections, e.g.  {{1,2},{1,4},{2,3},...} */
	  private double[][] BranchIndex;            
	  														/** (KxK) Negative Susceptance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
	  private double[][] NSM;     
	  														/** (KxK) Bus Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */ 
	  private double[][] BAM;       
	  														/** Reduced Bus Admittance Matrix <br> @see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */
	  private double[][] RBAM;
	  														/**(NxK) Adjacency Matrix <br>@see DC-OPF.IEEEPES2007.JSLT.pdf (pg 19)  */
	  private double[][] AM;         
	  														/** Temp matrix to form DA matrix*/
	  private double[][] dam;  
	  														/** Temp matrix to form RA matrix */
	  private double[][] ram;       
	  

	  														/** A factory producing dense matrices 1D */
	  private DoubleFactory1D Matrix1d = DoubleFactory1D.dense; 
	  														/**A factory producing dense matrices 2D*/
	  private DoubleFactory2D Matrix2d = DoubleFactory2D.dense; 
	  														/** Func is for using Colt's methods e.g. diagonal(), identity(), etc. */
	  private cern.jet.math.Functions Func = cern.jet.math.Functions.functions;


//C-O-N-S-T-R-U-C-T-O-R------------------------------------------------------------------------------------------------------------------------------  
//public OptimalPowerFlow(double[][] nodeData, double[][] branchData, double [] loadprofile, double[][] tsupplyoffer){
	  public OptimalPowerFlow(){
	  /*	
	    NodeDataA =     nodeData;
	    BranchDataA =   branchData;
	    LoadProfileA =  loadprofile;
	    TSupplyOfferA = tsupplyoffer;
*/	    
	   TSupplyOfferA = new double[5][4];
	   
			TSupplyOfferA[0][0] = 14.0;
			TSupplyOfferA[0][1] = 0.005;
			TSupplyOfferA[0][2] = 0.0;
			TSupplyOfferA[0][3] = 110.0;
			TSupplyOfferA[1][0] = 15.0;
			TSupplyOfferA[1][1] = 0.006;
			TSupplyOfferA[1][2] = 0.0;
			TSupplyOfferA[1][3] = 100.0;	
			TSupplyOfferA[2][0] = 25.0;
			TSupplyOfferA[2][1] = 0.01;
			TSupplyOfferA[2][2] = 0.0;
			TSupplyOfferA[2][3] = 520.0;
			TSupplyOfferA[3][0] = 30.0;
			TSupplyOfferA[3][1] = 0.012;
			TSupplyOfferA[3][2] = 0.0;
			TSupplyOfferA[3][3] = 200.0;
			TSupplyOfferA[4][0] = 10.0;
			TSupplyOfferA[4][1] = 0.007;
			TSupplyOfferA[4][2] = 0.0;
			TSupplyOfferA[4][3] = 600.0;	
			
		LoadProfileA = new double[3];
		
			LoadProfileA[0] = 350;
			LoadProfileA[1] = 300;
			LoadProfileA[2] = 250;

	    NodeDataA = new double[1][2]; 
	    
	    	NodeDataA[0][0] = 5.0;
	    	NodeDataA[0][1] = 0.05;
	    	
	    BranchDataA	= new double[6][4];
	    
	    	BranchDataA[0][0] = 1;
	    	BranchDataA[0][1] = 2;
	    	BranchDataA[0][2] = 250;
	    	BranchDataA[0][3] = 0.0281;
	    	BranchDataA[1][0] = 1;	
	    	BranchDataA[1][1] = 4;
	    	BranchDataA[1][2] = 150;
	    	BranchDataA[1][3] = 0.0304;
	    	BranchDataA[2][0] = 1;
	    	BranchDataA[2][1] = 5;
	    	BranchDataA[2][2] = 400;
	    	BranchDataA[2][3] = 0.0064;
	    	BranchDataA[3][0] = 2;
	    	BranchDataA[3][1] = 3;
	    	BranchDataA[3][2] = 350;
	    	BranchDataA[3][3] = 0.0108;
	    	BranchDataA[4][0] = 3;
	    	BranchDataA[4][1] = 4;
	    	BranchDataA[4][2] = 240;
	    	BranchDataA[4][3] = 0.0297;
	    	BranchDataA[5][0] = 4;
	    	BranchDataA[5][1] = 5;
	    	BranchDataA[5][2] = 240;
	    	BranchDataA[5][3] = 0.0297;
	    	
	        BranchDataA2	= BranchDataA;
	    	
	    nNodes = (int)NodeDataA[0][0];
	    nBranches = BranchDataA.length;
	    	
	    BranchIndex = new double[nBranches][2];
	    	
	  	 	K =  nNodes;          
	  	 	N =  nBranches;           
	  	 	I =  TSupplyOfferA.length;      
	  	 	J =  LoadProfileA.length;
	  	 	
	  	 GenAtNodeA = new int[I];
	  	 	
	  	 	GenAtNodeA[0] = 1 ;
	  	 	GenAtNodeA[1] = 1 ;
	  	 	GenAtNodeA[2] = 3 ;
	  	 	GenAtNodeA[3] = 4 ;
	  	 	GenAtNodeA[4] = 5 ;
	    	
	  	 LSEAtNodeA = new int[J];
	  	  	LSEAtNodeA[0] = 2;
	  	  	LSEAtNodeA[1] = 3;
	  	  	LSEAtNodeA[2] = 4;
	  	  	
	    InitParameters init = new InitParameters();
	    
	    
    for(int n=0; n<BranchDataA.length; n++){
        // Convert Branch Capacity from SI to PU
    	BranchDataA[n][2] = BranchDataA[n][2]/InitParameters.getBaseS();
    	
        // Convert Reactance from SI to PU, x(pu) = x/Zo = x/(Vo^2/So) = (x*So)/Vo^2
    	BranchDataA[n][3] = (BranchDataA[n][3]*InitParameters.getBaseS())/(InitParameters.getBaseV()*InitParameters.getBaseV());
      }
        //BranchDataA = Support.correctRoundingError(BranchDataA);
	
// SI to PU conversion for supply offer and load profile
      double sumCapL=0.0;
      double sumCapU=0.0;
      
      // Convert hourly Load Profile from SI to PU
      for(int n=0; n<LoadProfileA.length; n++){
      LoadProfileA[n] = LoadProfileA[n]/InitParameters.getBaseS();}
      //  LoadProfileA = Support.correctRoundingError(LoadProfileA);
      
      
      for(int i=0; i<TSupplyOfferA.length; i++){
       sumCapL+=TSupplyOfferA[i][CapLInd];
       sumCapU+=TSupplyOfferA[i][CapUInd];
       //System.out.println("\n supplyOfferByGen-> dMaxGenCapacity "+supplyOfferByGen[i][CAP_MAX]+" for GenCo :"+i);

       // Convert A from SI to PU-adjusted
       TSupplyOfferA[i][aInd] = TSupplyOfferA[i][aInd]*InitParameters.getBaseS();

       // Convert B from SI to PU-adjusted
       TSupplyOfferA[i][bInd] = TSupplyOfferA[i][bInd]*InitParameters.getBaseS()*InitParameters.getBaseS();

       // Convert CapMin from SI to PU
       TSupplyOfferA[i][CapLInd] = TSupplyOfferA[i][CapLInd]/InitParameters.getBaseS();

       // Convert CapMax from SI to PU
       TSupplyOfferA[i][CapUInd] = TSupplyOfferA[i][CapUInd]/InitParameters.getBaseS();
          
     }
     // TSupplyOfferArray = Support.correctRoundingError(TSupplyOfferArray);

	
	  TSupplyOfferMatrix = new DenseDoubleMatrix2D(TSupplyOfferA);
	  BranchIndex = new double[nBranches][2];  //e.g., {{1,2},{1,4},{2,3},...}
	
      NodeDataM = new DenseDoubleMatrix2D(NodeDataA);
      BranchDataM = new DenseDoubleMatrix2D(BranchDataA);
      
      
		
	  solveOPF( );
		
	}
	
//------------------------------------------------------------------------------------------------------------------------------------------    	
	private void solveOPF(){
		
		
	    setBranchIndex();
	    
	    formGM();
	    formaM();
	    formCeqM();
	    formbeqM();
	    formCiqM();
	    formbiqM();
		
	    qpj = new QuadProgJ(GM,aM,CeqM,beqM,CiqM,biqM);
	    boolean bHaveSolution=qpj.getIsFeasibleAndOptimal();

	    Commitment 	= new double[I];   // in MWs
	    VoltAngleR 	= new double[K-1];  // in radians
	    VoltAngleD  = new double[K-1]; // in degress    
	    IneqMultiplier = new double[2*N+2*I];
	    LMP         = new double[K];
	    minTVC = 0;  // in $/h
	    BranchFlow = new double[N]; //in MWs
	    double[] FullVoltAngle = new double[K]; // including delta_1 = 0
	    sumSquaredAngleDifference = 0;
	    
	    if(bHaveSolution) { // QuadProgJ has a solution
	        // DC-OPF solution for (p_{G1},...,p_{GI}) in SI
	        for(int i=0; i<I; i++){
	        	Commitment[i] = qpj.getMinX()[i]*InitParameters.getBaseS();
	        	
	        	System.out.println("Commitment" + Commitment[i]);
	        	
	        }

	        // DC-OPF solution for (delta_2,...,delta_K)
	        for(int k=I; k<I+K-1; k++){
	            VoltAngleR[k-I] = qpj.getMinX()[k];  // voltAngle in radians
	         }   

	        // Convert voltage angle from radian to degree
	        for(int k=1; k<K-1; k++){
	          VoltAngleD[k] = (VoltAngleR[k]*180)/Math.PI; // volt angle in degrees
	          
	          System.out.println("V Angle" + VoltAngleD[k]);
	          
	        }

	        for(int j=0; j<2*N+2*I; j++){
		      IneqMultiplier[j] = qpj.getIneqMultipiers()[j]/InitParameters.getBaseS();
		      
		      System.out.println("IneqMultiplier" +  IneqMultiplier[j]);
		      
		    }
	        
	        



	        for(int i=0; i<I; i++){
	          minTVC = minTVC + (A.get(i)/InitParameters.getBaseS())*Commitment[i]
	          +(B.get(i)/(InitParameters.getBaseS()*InitParameters.getBaseS()))*Commitment[i]*Commitment[i];
	        
	          System.out.println("TVC" +  minTVC);
	        
	        }

	        for(int k=1; k<K; k++){
	          FullVoltAngle[k] = VoltAngleR[k-1];
	        }
	        for(int n=0; n<N; n++){
	            BranchFlow[n] = (1/BranchDataA2[n][3])*(FullVoltAngle[(int)BranchIndex[n][0]-1]
	                         - FullVoltAngle[(int)BranchIndex[n][1]-1])*InitParameters.getBaseS();

	          sumSquaredAngleDifference = sumSquaredAngleDifference
	              + Math.pow((FullVoltAngle[(int)BranchIndex[n][0]-1]
	                          - FullVoltAngle[(int)BranchIndex[n][1]-1]),2);
	    
	        
	          System.out.println("BranchFlow" +  BranchFlow[n]);
	        
	        }
	        
	        // lmp: locational marginal prices in SI
	        for(int k=0; k<K; k++){
	          LMP[k] = qpj.getEqMultipliers()[k]/InitParameters.getBaseS();
	        
	          System.out.println("LMP" +  LMP[k] );
	        
	        }
        }
	}
	
//------------------------------------------------------------------------------------------------------------------------------------------    	
/** Forms Voltage Angle Difference Weight Matrix
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */	 
	
	private void formGM() {
		
	    formVADWM();
	    formRVADWM();
		 
		 B = new DenseDoubleMatrix1D(TSupplyOfferMatrix.viewColumn(bInd).toArray()); // Copies converted (/pu) values of b from supply offer into array B
		 UM = new DenseDoubleMatrix2D(Matrix2d.diagonal(B.assign(Func.mult(2))).toArray()); // Matrix U has 2*B diagonal and 0 for other elements
		 WrrM = new DenseDoubleMatrix2D(RVADWM);
		 GM = new DenseDoubleMatrix2D(I+K-1,I+K-1);
		 GM.assign(Matrix2d.composeDiagonal(UM,WrrM));
	 }
	 
	 
	 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Voltage Angle Difference Weight Matrix
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 16) */ 
	 private void formVADWM(){
		 
		 VADWM = new double[nNodes][nNodes];
		 
		    for(int n=0; n<nBranches; n++){
		      VADWM[(int)BranchIndex[n][0]-1][(int)BranchIndex[n][1]-1]
		          = -2*NodeDataA[0][PenCoefInd];  //NOTE: there should be a factor 2 in front of penaltyCoeff
		      VADWM[(int)BranchIndex[n][1]-1][(int)BranchIndex[n][0]-1]
		          = -2*NodeDataA[0][PenCoefInd];
		    }
		    for(int i=0; i<nNodes; i++){
		      for(int j=0; j<nNodes; j++){
		        if(j==i){
		          for(int k=0; k<nNodes; k++){
		            if(k!=i){
		            	VADWM[i][j] = VADWM[i][j] - VADWM[i][k];
		            }
		          }
		        }
		      }
		    }
		 //   VADWMatrix = Support.correctRoundingError(VADWMatrix); // Need a support Class
		  }
	 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Reduced Voltage Angle Difference Weight Matrix by deleting of VADWMatrix first row and first column
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 17) */ 	
   private void formRVADWM(){
	   
	   RVADWM = new double[nNodes-1][nNodes-1];
	   
		    for(int i=0; i<nNodes-1; i++){
		      for(int j=0; j<nNodes-1; j++){
		    	  RVADWM[i][j] = VADWM[i+1][j+1];
		      }
		    }
		  }
   
   
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Matrix of a coefficients from the supply offer
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) 
*  alpha*P+beta*P^2
*  */ 
   private void formaM(){
     A = new DenseDoubleMatrix1D(TSupplyOfferMatrix.viewColumn(aInd).toArray());

     aM = new DenseDoubleMatrix1D(I+K-1);
     aM.viewPart(0,I).assign(A); // viewPart(fromIndex, width)       // ?????????????????????????????
   }
   
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms  Adjacency Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */ 
      private void formAM(){  
    	  
    	  AM = new double[nBranches][nNodes];
    	  
    	    for(int n=0; n<nBranches; n++){
    	      for(int k=0; k<nNodes; k++){
    	        if(k == BranchIndex[n][0]-1){
    	        	AM[n][k] = 1;
    	        }
    	        else if (k == BranchIndex[n][1]-1){
    	        	AM[n][k] = -1;
    	        }
    	        else{
    	        	AM[n][k] = 0;
    	        }
    	      }
    	    }
    	  }
      
      
      
//------------------------------------------------------------------------------------------------------------------------------------------          
/** Forms Reduced  Adjacency Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */  
      private void formRAM(){
    	  
    	  ram = new double[nBranches][nNodes-1];
    	  
    	    DoubleMatrix2D amx = new DenseDoubleMatrix2D(AM);
    	    ram= amx.viewPart(0,1,nBranches,nNodes-1).toArray();
    	    RAM = new DenseDoubleMatrix2D(ram);
    	  }
      
 
//------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms II Matrix which represents K Nodes(rows)  and I Generators(columns) located on them.
 * <f/> If generator Ii is located at node Ki then the matrix element has value of 1 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */    
   private void formIIM(){
   IIM = new DenseDoubleMatrix2D(K,I);
   for(int k=0; k<K; k++){
     for(int i=0; i<I; i++){
       if(GenAtNodeA[i]==k+1){
         IIM.set(k,i,1);  // change set to setQuick() later
       }}}}

//------------------------------------------------------------------------------------------------------------------------------------------      
/** Forms Bus Admittance Matrix.
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */    
    private void formBAM(){
    	  
    	BAM = new double[nNodes][nNodes];
    	
    	  for(int i=0; i<nNodes; i++){
    	      for(int j=0; j<nNodes; j++){
    	        if(j==i){
    	          for(int k=0; k<nNodes; k++){
    	            if(k!=i){
    	            	BAM[i][j] = BAM[i][j] + NSM[i][k];
    	            }
    	          }
    	        }
    	        else
    	        	BAM[i][j] = - NSM[i][j];
    	      }
    	    }
    	  
    	  
      }

//------------------------------------------------------------------------------------------------------------------------------------------       
/** Forms Reduced Bus Admittance Matrix.
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
    private void formRBAM(){
    	
    	RBAM = new double[nNodes-1][nNodes];
    	
    	    DoubleMatrix2D na = new DenseDoubleMatrix2D(BAM);
    	    RBAM = na.viewPart(1,0,nNodes-1,nNodes).toArray();
    	    }
   
//------------------------------------------------------------------------------------------------------------------------------------------        	    
 /** Forms Negative Susceptance Matrix.
  *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 18) */
    private void formNSM(){
    	
    	NSM = new double[nNodes][nNodes];
    	
    	for(int n=0; n<nBranches; n++){
    		NSM[(int)BranchIndex[n][0]-1][(int)BranchIndex[n][1]-1]
    		          = 1/BranchDataM.viewColumn(ReactInd).toArray()[n];
    		NSM[(int)BranchIndex[n][1]-1][(int)BranchIndex[n][0]-1]
    		          = 1/BranchDataM.viewColumn(ReactInd).toArray()[n];
    		    }
    		  }
    	    
 //------------------------------------------------------------------------------------------------------------------------------------------    
/** Forms Equality Constraint Matrix
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 19) */ 
   // CeqTranspose = (II, -Br'); Ceq = CeqTranspose'; where Br' is rBusAdm here
   private void formCeqM(){
	   
	    formNSM();
	    formBAM();
	    formRBAM();
	    formIIM();
       
	   DoubleMatrix2D rBAM = new DenseDoubleMatrix2D(RBAM);
       DoubleMatrix2D[][] parts = {{ IIM, rBAM.viewDice().assign(Func.neg)}};
       DoubleMatrix2D CeqTMatrix = new DenseDoubleMatrix2D(K,I+K-1);
       CeqTMatrix.assign(Matrix2d.compose(parts));
       CeqM = new DenseDoubleMatrix2D(I+K-1,K);
       CeqM.assign(CeqTMatrix.viewDice());
   
   }   
   
   
//------------------------------------------------------------------------------------------------------------------------------------------    

// FDemand = someFunction(atNodeByLSE, loadProfile); beq = FDemand
/** Forms an Equality Constraint Vector 
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */
       private void formbeqM(){

    	   beqM = new DenseDoubleMatrix1D(K);
         for(int k=0; k<K; k++){
           double L=0;
           for(int j=0; j<J; j++){
             if(LSEAtNodeA[j]==k+1){
               L = L + LoadProfileA[j];
             }
           }
           beqM.set(k,L);
         }
     
   }   

//------------------------------------------------------------------------------------------------------------------------------------------	 
       // Ciq matrix formulation
       //
       // (1) FIXED DEMAND CASE:
       // CiqTranspose = {{Oni, Z*rAdj},{Oni, -Z*rAdj},{Iii, Oik},{-Iii, Oik}};
       // Ciq = CiqTranspose'
    
/** Forms an Inequality Constraint Matrix Ciq 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 20) */     

       private void formCiqM(){

    	   formOtM();
   	       formOpM();
   	       formDAM();
   	       formAM();
   	       formRAM();
   	       formIpM();
    	   
    	   
       
           DoubleMatrix2D[][] parts = {
           { OtM,   DAM.copy().assign(Func.neg).zMult(RAM,null)},  // (1) changed upside-down 	   
           { OtM,                  DAM.zMult(RAM,null)         },  // (2) regarding AMES       
           { IpM,                             OpM                   },
           { IpM.copy().assign(Func.neg),     OpM                  }};
           
           
           DoubleMatrix2D CiqTMatrix = new DenseDoubleMatrix2D(2*N+2*I,I+K-1);
           CiqTMatrix.assign(Matrix2d.compose(parts));
           CiqM = new DenseDoubleMatrix2D(I+K-1,2*N+2*I);
           CiqM.assign(CiqTMatrix.viewDice()); 
       }
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms Diagonal Admittance Matrix, where Bkm=1/Xkm 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */       
       private void formDAM(){
    	   
    	   dam = new double[nBranches][nBranches];
    	   
    	   dam = Matrix2d.diagonal(BranchDataM.copy().viewColumn(ReactInd).assign(Func.inv)).toArray();
    	   DAM = new DenseDoubleMatrix2D(dam);}
    	    //NOTE: Have to keep .copy(), otherwise reactance will be 1/reactance
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms O an N × I zero Matrix 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
       private void formOtM(){
    	   
    	   OtM   = new DenseDoubleMatrix2D(N,I);     	   
       }
 
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms O an (I x K-1) zero Matrix 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
        private void formOpM(){
           	   
        	OpM = new DenseDoubleMatrix2D(I,K-1); 
        }
       
//------------------------------------------------------------------------------------------------------------------------------------------	 
/** Forms I×I Identity Matrix 
 *  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
        private void formIpM(){
                   	   
        	IpM = new DenseDoubleMatrix2D(I,I).assign(Matrix2d.identity(I));
        }  
        
        
 //------------------------------------------------------------------------------------------------------------------------------------------	 
	 private void setBranchIndex(){
		    BranchIndex = BranchDataM.viewPart(0,FromInd,nBranches,2).toArray();
		    }
	 


//------------------------------------------------------------------------------------------------------------------------------------------
// biq = (-pU, -pU, capL, -capU)
/** Forms (2N+2I)×1 the associated inequality constraint vector 
*  @see <f/> DC-OPF.IEEEPES2007.JSLT.pdf (pg 21) */
	 private void formbiqM(){

	
	 
	 BCapM =new DenseDoubleMatrix1D(BranchDataM.viewColumn(BCapInd).toArray()); // !!! Check if it's working	 
	 CapTLM = new DenseDoubleMatrix1D(TSupplyOfferMatrix.viewColumn(CapLInd).toArray());
	 CapTUM = new DenseDoubleMatrix1D(TSupplyOfferMatrix.viewColumn(CapUInd).toArray());

	 
     DoubleMatrix1D[] parts = {BCapM.copy().assign(Func.neg), BCapM.copy().assign(Func.neg),
    		                  CapTLM, CapTUM.copy().assign(Func.neg)};
    
    
     biqM = new DenseDoubleMatrix1D(2*N+2*I);
     biqM.assign(Matrix1d.make(parts));
    
     }
	 
	 public void check(){
		 
		 System.out.println("Proverka na prochnost");
	 }
    }












